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Abstract 

We estimate the momentum diffusion coefficient of a heavy quark within a pure SU(3) plasma 
at a temperature of about l.hTc. Large-scale Monte Carlo simulations on a series of lattices 
extending up to 192^ x 48 permit us to carry out a continuum extrapolation of the so-called 
colour-electric imaginary-time correlator. The extrapolated correlator is analyzed with the 
help of theoretically motivated models for the corresponding spectral function. Evidence for a 
non-zero transport coefficient is found and, incorporating systematic uncertainties reflecting 
model assumptions, we obtain k = (1.8 — 3.4) T^. This implies that the “drag coefficient”, 
characterizing the time scale at which heavy quarks adjust to hydrodynamic flow, is = 
(1.8 — 3.4)(rc/r)^(M/1.5GeV) fm/c, where M is the heavy quark kinetic mass. The results 
apply to bottom and, with somewhat larger systematic uncertainties, to charm quarks. 


^In memoriam (27 October 1956 - 1 June 2015). 



1. Introduction 


Within linear response theory the rate at which a system relaxes towards local thermal 
equilibrium is characterized by quantities known as transport coefficients. Different transport 
coefficients parametrize different types of perturbations. If we focus on a conserved particle 
number, such as a quark flavourj^ which is initially distributed unevenly, such as in a broad 
jet cone, then the relevant transport coefficient is the diffusion coefficient. In QCD, there is 
a separate diffusion coefficient related to light flavours, and to heavy flavours such as charm 
and bottom. A closely related quantity is the electrical conductivity, which can be expressed 
as a weighted sum over the flavour diffusion coefficients. 

The determination of transport coefficients related to strong interactions at temperatures 
of a few hundred MeV is an important goal for lattice QCD. It is a challenging task, given that 
lattice QCD is formulated in a Euclidean spacetime whereas transport coefficients are real¬ 
time quantities (for a review, see ref. [T]). Nevertheless, motivated by the ongoing heavy ion 
collision program, large-scale efforts have been undertaken. For example, for the light-quark 
contribution to electrical conductivity, recent works can be found in refs. 12H7]. 

The focus of the present study is the diffusion coefficient associated with heavy quarks. 
It has been one of the major qualitative discoveries of the heavy ion collision program at 
RHIC and LHC that charm quarks appear to flow about as efficiently as light quarks do (see, 
e.g., refs. isHni] and references therein). That a flow develops is an example of a relaxation 
towards local thermal equilibrium. It is then a theoretical challenge to explain this behaviour 
from the laws of QCD m- A next-to-leading order (NLO) computation in perturbation the¬ 
ory indicates the presence of a large correction towards strong interactions [12], and strong 
interactions have also been observed in AA = 4 Super-Yang-Mills theory |13H16] . Furthermore 
classical lattice gauge theory simulations and analyses in the confined phase [18II22] are 
consistent with strong interactions, and various other approaches are being pursued in the 
same vein [23ti29] (for a review, see ref. [30]). Heavy quark diffusion also happens to pose 
an ideal ground for more general theoretical investigations of non-equilibrium thermodynam¬ 
ics [HI]. In any case, ultimately the problem needs to be addressed with lattice simulations. 

For M ^ ttT, where M is the heavy-quark “kinetic mass”, the lattice determination of the 
heavy quark diffusion coefficient can be reduced to a purely gluonic measurement [32]. Here 
we report the final results of a multiyear study of the relevant observable in the deconfined 
phase of pure SU(3) gauge theory. It has been demonstrated a while ago that, with advanced 
numerical methods, a signal can be obtained at a fixed lattice spacing [33II35] . However the 
issues of renormalization, taking the continuum limit, and analytic continuation had not been 
brought into conclusion. Even though further improvements are needed and can be foreseen 
on all of these fronts, the purpose of the current paper is to present an analysis which offers 

^Weak interactions play no role within the lifetime < 20 fm/c of a fireball generated in a heavy ion collision. 
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Table 1: The lattices included in the current analysis. Conversions to units of tg are based on ref. [36] 
for the improved Wilson discretization (“imp”); on refs. [36l|37] for the clover discretization (“clov”) 
of the observable defining tg [35]; and on ref. [33] for rg [33]. Note however that the clover case is not 
well represented by our ansatz [36]: y^/d.o.f. ~ 42. Conversions to Tc are based on ref. [36] . 



Figure 1: One possible discretization of eq. (12.1|) [32] . Different techniques have been used for 
improving on the statistical signal originating from the links denoted with thick lines, and from the 
electric field insertions delineated with the vertical lines (cf. the text). 

a “minimal” practical answer to the main open points. 

The plan of this paper is the following. After introducing the observable and reviewing 
the techniques that have been used for its determination and the data sets that have been 
collected (sec. Ej), we carry out a continuum extrapolation in sec. [3j The simple structure 
of the associated spectral function, as revealed by previous theoretical works, allows us to 
attack the difficult problem of analytic continuation through tightly constrained models, as 
well as a variant of the so-called Backus-Gilbert method (sec. 0]). Our final results and some 
future prospects are presented in sec. (5] 


2. Measurements 


Using Heavy Quark Effective Theory, the force felt by a heavy quark as it propagates through 
a gluon plasma can be related to a “colour-electric correlator” [32] (cf. also ref. [T5]l. 
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Figure 2: Left: Lattice data after perturbative renormalization and the use of tree-level improved 
distances. Right: Illustrations of continuum extrapolations at selected values of tT (our final errors 
are based on a jackknife analysis and differ somewhat from those shown here.) 


where gEi denotes the colour-electric field, T the temperature, and U{t 2 ',ti) a Wilson line in 
the Euclidean time direction. The discretization of this correlator is not unique; we employ 
the proposal of ref. [32], illustrated in fig. (U 

We have measured the discretized correlator in quenched lattice QCD, employing the stan¬ 
dard Wilson gauge action, at a temperature corresponding to about 1.5Tc (the parameters of 
the simulations are shown in tabled]). In order to obtain a signal for our observable, advanced 
statistical error reduction techniques are required. As has been discussed in more detail in 
refs. [MIHO] . we have employed 1000 additional multilevel updates [^TM2] for the electric field 
insertions, and link integration [431144] for the straight lines between them. Moreover, in order 
to reduce discretization effects, the imaginary-time separations are tree-level improved [33] 
in analogy with the procedure previously used in other contexts [3^142] . 

The resulting correlators are illustrated in fig. [2](left). For better visibility, the correlators 
have been normalized to (cf. ref. [32]) 




,(r) = 


cos^(7rTT) 

sin^(7rrT) 3sin^(7rrT) 


+ 


( 2 . 2 ) 


which diverges as l/(7r^r^) at r <C /3. The correlators have also been multiplied by a renor¬ 
malization factor as will be explained in sec. [31 As discussed in more detail in ref. [34] . 
volume dependence lies below statistical uncertainties in our measurements. In the following, 
we therefore consider a fixed spatial extent in units of the temperature, Lg = 4/r. 
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2.3001(44) 
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3.4005(161) 
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23/48 
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Table 2: The continuum-extrapolated correlator, normalized to eq. (12.21) . The errors are statistical 
and result from linear extrapolations in l/Al^ to the continuum limit (cf. the text). 


3. Continuum extrapolation 

Given the lattice data at tree-level improved distances, we carry out an extrapolation to 
the continuum limit. For this the lattice-measured correlator needs to be multiplied by a 
renormalization factor: 


. (3.1) 

A 1-loop perturbative computation yields [35] 

^E,pen = 1 + 0.079 X A + o(o.079 xA)', (3.2) 

where /3 q = Q/qq is the coupling of the plaquette term in the Wilson action (cf. tabled]). 
Ultimately renormalization should be carried out non-perturbatively, however the small coef¬ 
ficient of the 1-loop term in eq. (13.2p suggests that perturbative renormalization should yield 
a reasonable approximation. 

Measured results at 4 or 5 different values of (depending on tT, cf. fig. [2|) are interpo¬ 
lated to the values of tT that are shown in table [2j At fixed tT we extrapolate the correlator 
in to the continuum limit. The procedure is illustrated in fig. [2Kright) for selected values 
of tT. The resulting continuum limit is shown in fig. [HKleft), and the results are tabulated in 
table [2] The results and errors of continuum-extrapolated correlation functions were obtained 
from a combined jackknife analysis. A covariance matrix for the continuum correlator was 
also estimated within this analysis (however the errors are strongly correlated even at large 
r-separations; the covariance matrix has very small eigenvalues and oscillatory eigenvectors, 
and therefore its inverse is of limited practical use in fitting). 
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4. Modelling the spectral function 


Given the data for the imaginary-time correlator, the next task is to constrain the corre¬ 
sponding spectral function. The relation of a spectral function p^{uj) to the corresponding 
imaginary-time correlator G^(t) reads 


G,(r) = 


Jo ^ 


cosh[a;(^ — r)] 


(4.1) 


Even though an inversion of this relation is possible in principle [46] (after the subtraction 
of short-distance singularities m), the problem is ill-posed in practice: large variations of 
Pg may lead to small changes of G^. Therefore, it is important to constrain the allowed 
form of Pg from general considerations. Here we do this by fixing the functional form of 
Pg at small (w <C T) and large frequencies (cu S> T). Subsequently theoretically motivated 
interpolations between the two limits are proposed. We consider the spatial volume to be 
infinite and assume pg to be a smooth function of cu, as is generally the case in an interacting 
thermal system. 


4.1. IR and UV asymptotics 


In the infrared (IR) regime (w <C T), the heavy quark momentum diffusion coefficient can be 
defined as [32] 


K = lim 

cj —>-0 


2rpg(w) 


u 


(4.2) 


The approach to this limit appears to be smooth: resummed perturbative computations 
numerical simulations within classical lattice gauge theory m, as well as strong-coupling 
computations in analogous theories [TSlIlH] suggest that pg has no transport peak but is 
rather a monotonically increasing function. Therefore, we define the infrared asymptotics 
through the simplest form consistent with eq. dMj): 



(4.3) 


Consider then the ultraviolet (UV) regime (w ^ T). Thanks to asymptotic freedom, the 
UV behaviour of the spectral function can be computed in perturbation theory. Denoting by 
the QCD gauge coupling renormalized in the MS scheme, and = a^/vr = p^/(47r^), the 
result has the structure 




J r=o 


+ 0 




UJ 


(4.4) 


The leading thermal correction is consistent with the pattern expected from the Operator 
Product Expansion |l9] (for Vf = 0 the coefficient of this correction is negative [50]). The 
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vacuum part has the form 
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where Cp = (JV^ — l)/(2Nc) = 4/3 and £ = ln(/x^/a;^), with /i denoting the renormalization 
scale. The three first coefficients read |50j 
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(4.6) 

(4.7) 

(4.8) 


where Nf denotes the number of light dynamical quarks (iVf = 0 in our study). The coeffi¬ 
cient Tgg and higher-order terms are presently unknown. However, the general structure of 
eq. (14.5p . together with the knowledge that Pp requires no renormalization in dimensional 
regularization [32], is sufficient for determining the asymptotics of Pp. 

Indeed, suppose that we choose the renormalization scale as /i = w for cu ^ Then 

£ = 0, ~ ln“^(w/Ai^), and we get 


/ \ i{o,) / \ 

^ „ = '?>uv(w) 


J T=0 


1 + 0 


1 


ln(a;/AMs)y_ 


(4.9) 


where we have defined 

A<P), \ - a^iPjCpUJ^ _ _ I rT.^ 

</>uv(^^) = -^= max(w,vrT) . (4.10) 

OTT 

Formally, we can reduce the correction to be of quadratic rather than linear order in ln“^(a;/Aj^) 
by including r 2 Q in the asymptotics, so let 


^>uv(w) = </’uv(^^) 


1 + {r20 + r2i£) a,{pj 


(4.11) 


However, since the convergence of an expansion proceeding in inverse logarithms could be 
slow, and since our knowledge of a;/Aj^ = (w/T) X (r/A]y[g) is imperfect due to uncertainties 
in T/Aj^ [Ml, we treat eqs. (I4.10p and (|4.11l) on an equal footing in the following. For 
evaluating and ag{p^), we have used 4-loop running [5T] . 


4.2. Interpolations 

Equations (|4.3I1 and (|4.10l) determine the limiting behaviours of the spectral function. For 
a; <C T, given the flatness observed in non-perturbative simulations m, we expect corrections 
to to be given by a convergent power series in lo. For uj A> T the corrections are suppressed 
by inverse logarithms of oj/T, and it is important to account for these corrections. In order 
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to incorporate the different types of corrections in the two regimes, we map the interval 
Ld G (0, oo) to the interval (0,1) by introducing 


( UJ \ X 

l + ^)€(0,oo), y=—— g(0, 1). 

ttT / 1 + X 


For small y we have 


OJ 


— , uj<^T 
TTi 


whereas around the other end we get 


1-y 


1 


X ln[a;/(7rT)] 




(4.12) 


(4.13) 


(4.14) 


So, we can capture corrections to the asymptotics by a function in y which vanishes at y = 0 
and y = 1 and is a polynomial or a power series in between. A convenient choice is to employ 
trigonometric functions, 

4“^ (y) = sin(7rny) . (4.15) 

Another possibility, in principle preferable if we know that corrections are only quadratic in 
the variables in eqs. (14.1311 . (14.141) (cf. the discussion preceding eq. (|4.1ip i. is 

e^^\y) = sin(7ry) sin(7rny) . (4-16) 


With these bases we are led to define general models (// G {a, /3}, i G {a, b}), 

^max 

model 1 : = l+^c„e2 "^(y) (/>ij^(w) + 4v(^^) : 

n=l 

“^max j - 

model 2 : ^ c„ei^)(y) y 4 ir(w)]^ + [4 *4^^)]' 


n=l 


(4.17) 

(4.18) 


The difference between these interpolations is that the latter imposes a more rapid crossover 

from the IR to the UV asymptotics, for instance for large u) we have (t’m + ~ 4’vv + 

4^ ~ lo^+T^/co in qualitative accordance with the functional form expected from eq. (14.41) rl 
^<Puv 

Finally, we also consider a simple 2-parameter ansatz separating the IR and UV regimes 
completely: 


model 3 : 


(3i)^ N — 
p^ '[u) = max 


</'ir(‘^)>c4v(^) 


(4.19) 


where c is treated as a free parameter, reflecting uncertainties in the renormalization factor 
in eq. (13.2|) (in practice we find c = 1.05(1) from the fit for i = a). 

®Note that, for Umax —t oo, such a construction can be viewed as a general parametrization of an arbitrary 
spectral function. We call our parametrizations “models” because the problems discussed in sec. l4.3l necessitate 
keeping Wmax relatively small and stabilizing the subsequent fits through additional input. 
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4.3. Fitting strategy 


The inversion of eq. m represents an ill-posed problem, so it should not come as a surprise 
that it is, in general, not possible to find a “stable” fit describing the data. In other words, the 
X^-function may possess an extremely shallow minimum, or multiple minima. Empirically we 
find, however, that these ambiguities are largely related to the UV behaviour of the spectral 
function. With some further input, the UV behaviour can be stabilized, yet the IR coefficient 
K in which we are mostly interested, remains fairly unaffected. 

Concretely, we have implemented two strategies for the fitting, which incorporate an im¬ 
plicit or explicit stabilization of the UV contribution. After defining 

where denote the errors from table [2l the following alternatives are considered: 

(i) In the first strategy, we search for a minimum of for a result based on eq. M, 
taking a jackknife sample as the measured result. The search is performed using the 
algorithm described in ref. [S2], taking k/T^ = 1.0, Cn = 0.0 as a starting point. The 
search is stopped after 200 iterations, by which time an excellent representation of the 
data has in general been reached (with y^/d.o.f. ~ 0.3 — 0.5 for a typical sample from 
the jackknife ensemble), with coefficients |c„| <0.2. 

(ii) In the second strategy, we stabilize the fit by imposing the constraint 

, (4-21) 

where we choose = 1000 T. This constraint imposes a relation between k and 

c„, so that there are only free parameters. We find that with this constraint, a 
unique minimum of can always be found, with y^/d.o.f. ~ 0.2 — 0.6. The practical 
search was performed with the routine described in ref. m- Because of the artificial 
constraint at large uj/T, the fit is carried out only to distances tT > 6/48. 

In the following we describe results from both procedures, demonstrating that they yield 
similar results for k, particularly for our preferred “model 2” (cf. eq. (I4.18p i. 

In the case of the 2-parameter “model 3” dehned in eq. (I4.19p . a global minimum of y^ 
is readily found without further input. The price to pay for this stability is that the model 
does not describe the data particularly well, having y^/d.o.f. « 17.5 for (j)^l (we do not show 
results for for which y^/d.o.f. « 68). 

“^We have also carried out tests with the covariance matrix, but in its full form this does not yield sensible 
results as alluded to above. If a sufficient infrared cutoff is imposed on the smallest eigenvalues, the results 
are consistent within errors with the procedure described here. 


^meas(R) ^model(u) 


• z 


(4.20) 










In addition to these fits, we have also made use of a variant of the Backus-Gilbert method 
(BGM) [MIES] (see refs. [561EZ! for previous applications in lattice QCD). The goal of this 
method is not to reconstruct the spectral function itself, but rather an averaged version 
thereof. With very precise data, the averaging kernel could be made optimally narrow in 
a certain sense. In practice, the finite precision of the data necessitates a regularization 
of an ill-defined matrix inversion; this is characterized by a parameter A for which we use 
A = rather than the theoretically optimal A = 1 (cf. ref. [57]). As a consequence, our 
estimate of the infrared limit of the spectral function amounts to a weighted average over the 
range 0 < cj < lOT. Fortunately, if there is little structure in this range, the result should 
be reasonable. Moreover it is possible to rescale the kernel in order to further remove known 
structures; we insert Py^{oj) = [p-^{ijj) / (j){ijj)\ 4>{oj) in eq. (j4.ip with 4>i^) = / tanh^(a;/3/4), 

and subsequently include (p in the kernel. To keep the matrix size manageable, only t/13 = 
(4 — 21)/48 were used for this analysis. If the final result for the averaged p^{oj) is inserted 
back into eq. (BU), it yields a representation of the original data with x^/d.o.f. ~ 41. 

We end by remarking that apart from these novel approaches, we have also applied the 
standard Maximum Entropy Method (MEM) to our problem. MEM requires the specification 
of a default model as input. If we take as the default model our model 1 or 2, which yield 
a Euclidean correlator agreeing with lattice data everywhere (within statistical errors), then 
MEM reproduces the default model as its output (within statistical errors). In other words, 
MEM does not help to constrain the result beyond our analysis. 


4.4. Estimation of k 

We now turn to estimating k. In sec. 14.21 we have introduced a large class of fit functions: 
8 separate models p^E^^\ ^ £ {1)2}, p G {q;,/3}, i G {a,^, each of which depends on a 
parameter In addition we have introduced two fitting strategies (cf. sec. 14 . 3p . In the 

following, we choose = 4 , 5 and demonstrate that, for a given model, both strategies 

yield similar results within error bars. (We have also used = 3 and = 6; these yield 
nothing qualitatively new.) The agreement is more significant for our preferred model 2 . The 
spread between the models is interpreted as an indication of the systematic uncertainty of 
our determination. Model 3a and the BGM approach are included as further crosschecks; in 
the BGM case the spread of results originating from different versions of covariance matrices 
(cf. caption of fig. HP is included in the error shown. 

The fits are carried out to samples from a jackknife ensemble, and the ensemble is used 
for determining statistical errors. Results from the different models are illustrated in figs. [3l 
ll]and[5](the fits in figs. [3] and H] are based on strategy (i)). Given that model 2 yields more 
stable results within the other variations, and that it is theoretically better justified than 
model 1 (cf. the discussion below eq. (I4.18D ). we make use of it in the following. Based on 
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Figure 3: Left: Continuum-extrapolated lattice data (cf. tabled]) and examples of model results, as 
described in the text. The results have been normalized to eq. (12.211 . Right: The differences between 
the data and the models. For better visibility, the data points corresponding to models laa and 3a 
have been displaced slightly. 


the central values from fig. [5] for models 2, 3a and BGM we estimate 


a/T^ = 1. 


3.4 


(4.22) 


This result is illustrated in fig. [5| with a grey band. Given that the result is dominated by 
systematic uncertainties which may be asymmetric but are essentially impossible to estimate, 
we refrain from citing a central value here@ 


5. Conclusions 

The purpose of this paper has been to estimate the heavy-quark momentum diffusion coef¬ 
ficient, defined through eqs. (EH), (SID and S2D- Compared with previous works |33H35j . 
we have carried out a continuum extrapolation of the imaginary-time correlator (cf. sec. [3D 
and discussed the systematics related to estimating the corresponding spectral function (cf. 
sec. HD- The final result of our analysis is given in eq. (|4.22p . It is remarkable that, despite 
the ill-posed nature of analytic continuation, our novel approach permits us to obtain rel- 

®It is remarkable that our result is consistent with previous ones obtained at a single lattice spacing and 
with very rough modelling of the spectral shape |341I35| . However, the uncertainties related to the continuum 
extrapolation and to analytic continuation are now much closer to being under control. 
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Figure 4: Examples of model spectral functions, compared with from eq. (14.101) . In the BGM 
case the curve shown has been obtained by replacing the full covariance matrix with a diagonal one 
possessing the errors shown in table [5J but we have obtained BGM results with full and IR-regularized 
covariance matrices as well; they are included in the average shown in fig. [5] 

atively stable results, and to strongly constrain the order of magnitude of the heavy quark 
momentum diffusion coefficient in the continuum limit. 

In the non-relativistic limit (i.e. for a heavy quark mass M ^ ttT) k is related to the 
diffusion coefficient D as D = 2T‘^j k, and to the drag coefficient as = k/{2MT), where 
M is a heavy quark kinetic mass [ 32 ] • The drag coefficient can also be interpreted as the 
kinetic equilibration time scale associated with heavy quarks: For a conversion 

to physical units, we use r^Tc = 0.7457(45) [36| and Tq = 0.47(1) fm [58|. With these 
conversions, our result k/T^ = 1.8 — 3.4 yields an estimate for the time scale associated with 
the kinetic equilibration of heavy quarks, 

r„.. = i = (1.8,,,3,4)(|) (5,1) 

Close to Tc, charm quark kinetic equilibration appears therefore to be almost as fast as that 
of light partons, for which a time scale ~ 1 fm/c is generally considered appropriate. For the 
diffusion coefficient we obtain 

nr = 0.59... 1.1 . (5.2) 

This can be compared with the values DT > 0.13 obtained for light quarks in quenched QCD 
in the continuum limit W- It would be “natural” for the D of heavy quarks to be of the 
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T~ 1.5T 



Figure 5: Fit results based on different models (cf. sec. l4.2l) and different fitting strategies (cf. sec. l4.3l) . 
In each case, the lower data point corresponds to = 4, the higher to = 5. For model 3a 
(cf. eq. (I4.19p l and BGM (cf. sec. 14.3p systematic errors are larger than those shown. The grey band 
illustrates our final estimate, given in eq. (I4.22|) and based on models 2, 3a and BGM. 

same order of magnitude but somewhat larger than that of light quarks, given that heavy 
quarks should feel slightly weaker interactions. 

It is also interesting to compare our result for k with an NLO perturbative computation |12j . 
That result is of the form k/T^ = {c^ + 020 ^^), where C;^,C 2 are coefficients given in 
ref. |12j [ci involves a logarithmic dependence on Og). In the absence of corrections of 
relative order as, it is not possible to estimate the renormalization scale at which as should 
be evaluated. Nevertheless, the result shown in fig. 3 of ref. [I2] agrees with our eq. (j4.22p if 
we set as = 0.20 — 0.26, which is in full accordance with the range generally used in heavy 
ion collision phenomenology. 

Many possible directions can be envisaged for future investigations. Improved statistical 
precision is crucial for moving towards model-independent analytic continuation m- Other 
temperatures than just T ~ l.STc should be considered. The determination of the renormal¬ 
ization factor in eq. (13.ip should be promoted to the non-perturbative level. It would be 
important to understand whether the heavy-mass limit is justified for charm quarks (or only 
for bottom quarks); this can in principle be studied by using the full relativistic formulation 
for measuring current-current correlation functions [59ll60j , even though then the structure of 
the spectral function is more complicated and analytic continuation is even more difficult to 
get under reasonable control. Finally, estimating effects from dynamical quarks is important 
for phenomenological applications. 
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